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I .  Introduction 


A  workshop  under  the  direction  of  Norbert  Peters  (Instituet  fuer  Allgeaeine 
Mechanik  RWTH  Aachen)  has  been  called  to  focus  on  numerical  Methods  that  have 
been  developed  to  solve  premixed  flame  problems  with  complex  chemistry  and 
transport  properties.  The  purpose  of  this  workshop  is  to  seek  to: 

(1)  Establish  the  difficulties  that  result  from  non-equal  diffusivities 
of  heat  and  matter  and  the  consequences  that  these  have  upon  the  accuracy  of 
the  solution; 

(2)  Compare  different  numerical  schemes  for  two  test  problems;  and 

(3)  Exchange  ideas  about  work  in  progress.  This  workshop  has  been  subti¬ 
tled  "The  Flat  Flame  Olympics,”  and  was  held  12-14  October,  1981  at  the 
University  of  Aachen.  Proceedings  of  the  entire  workshop  are  to  be  published. 1® 

This  report  covers  the  second  test  problem:  the  solution  of  the  steady, 
stoichiometric  hydrogen-air  flame  with  complex  chemistry.  Me  have  named 
it  the  Peters  B  problem.  (The  first  problem,  Peters  A,  is  detailed  in  another 
reportlb) .  The  kinetic  scheme  and  boundary  conditions  at  the  cooled  end  are 
specified  and  the  governing  equations,  boundary  conditions  at  the  burned  end, 
the  initial  conditions,  transport  properties  and  heat  capacities,  and  numerical 
scheme  were  left  to  the  choice  of  the  authors. 

This  report  then  documents  the  choices  made  and  presents  material  requested 
for  the  workshop  as  well  as  other  interesting  findings. 

II.  The  Flame  Equations 

The  derivation  of  the  conservation  equations  for  a  multicomponent  reacting 
ideal  gas  mixture  can  be  found  in  several  textslc_3.  Those  equations  that 
adequately  describe  a  one-dimensional,  laminar,  premixed  flame  that  propagates 
in  an  unbounded  medium  are: 

Overall  continuity 

pt  +  (pu}x  “  0 ;  (1) 


3 aNotes  on  Numerical  Fluid.  Mechanics,  Vol.  5,  N.  Peters  and  J.  Wamatz ,  Ed. , 

Vieweg- Vertag,  Pub. 

lbT.P.  Coffee ,  "Flat  Flame  Olympias:  Test  Problem  At"  BRL  Report  to  be  published. 

loHirschf elder ,  J.O.,  Curtiss ,  C.F.,  and  Bird ,  R.B. ,  Molecular  Theory  of  Oases 
and  Liquids ,  2nd  Ed.  (corrected  with  notes) ,  Wiley,  Few  Fork,  1964. 

2Bird,  R.B. ,  Stewart ,  W.S.  and  Light  foot,  E.S.,  Transport  Phenomena ,  Wiley, 

Hew  York,  1960. 

3 Williams,  F.A.,  Combustion  Theory,  Addison-Wesley ,  Reading,  MA,  1966. 


Continuity  of  species 


»<Vt  *  »“  tVx  ■  <»Wx  *  Vk 

and  conservation  of  energy. 


1,2,...,  N: 


ecpTt  +  PUCpTx  -  CXTx)x  -  j1  ♦  \ml  hk  (PYkVk): 


(2) 

C3) 


where  the  variables  are  defined  in  “.he  glossary.  The  effects  of  radiation, 
viscosity,  thermal  diffusion  and  body  forces  have  been  neglected  and  since  the 
burning  velocity  is  small  compared  to  the  local  speed  of  sound,  the  pressure  is 
taken  to  be  constant3-4.  Our  general  thrust  is  the  determination  of  species 
and  temperature  profiles  and  of  burning  velocities.  The  boundary  conditions 
are  the  following  (t>o) .  For  x  ■  -  «;  T  ®  T  and  Y.  =  Y.  ,  k=l,2,...,  N, 
and  for  x  *  +  ®  u  K  Ku 


Tx  *  (Vx  *  0  (k  a  1'2'  •••)• 


C4) 


In  order  to  avoid  solving  Eq  (1),  a  coordinate  f  is  introduced  such  that 


♦(x,t)  =  /  p(x  ,t)dx  . 


(5) 


Then  \ |»x  *p  and  i|»t  »  -pu  ♦  mQ(t), 


where  (PU)IX=Q  ”  m  (t)..  With  this  coordinate  change  equations  C2)  and  (3)  become 


and 


<Vt  +  W*  ‘  *  *kMk  » 

N 

=  cP"1{(pAV*  -  k'l  Wk 


(6) 


(7) 


respectively.  In  practice  convenient  dimensionless  forms  of  equations  (6)  and 
(7)  are  solved. 


The  burning  velocity,  S  ,  is  defined  as  the  velocity  of  the  flame  relative 
to  the  fluid  at  rest,  that  is,  at  "infinity."  The  value  of  Su  can  be  found 
from  the  steady-state  profiles  by  performing  a  coordinate  transformation  such 
that  the  flame  appears  stationary.  In  this  new  system  all  variables  are  inde¬ 
pendent  of  time,  specifically,  (Yfc)t  ■  ■  0 •  Then  from  Equation  (1)  we  have 

(pu)  *>  0  or  pu  ■  constant.  Take  any  one  of  the  equations  (2)  and  integrate 
overxthe  interval  (a,b)  to  obtain: 


pu  [Yk(b)  -  Yk(a)]  -  /Rb  RjMjto  -  pYkVk[J  , 

^  ,  x  '  a"  W*  -  >YkVk  \l  'a  P'1  W  -  pYkVk  >a 

"  pH  [Ykfl>)"-  "  p(--)  [Yk(b)  -  Yk(a))  ‘ 

4  Prietrom,  R.M.  and  Weetenberg ,  A. A.,  Combustion  Theory,  Addiaon-Wealey, 
Reading,  MA,  196 S. 


(8) 

(9) 
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.... 


In  this  coordinate  system  Su  ■  -u(-<°) • 

One  could  also  take  Bq.  (3)  and  integrate  it  over  (a,b)  to  find  an  expression 
similar  to  Eqs.  (8)  and  (?).  Specifically,  we  have 


u(-“) 


0  Jl  W  *  XTx  I." 


CIO) 


P(-)  /a  cpdT 

which  for  constant  c  reduces  to 
h  P  N 

V  kli  \w 


U(-“) 


x  +  XT. 


p(-)  cpiTb  -  Tar 


x  'a 


(ID 


For  the  transport  method  used  below  h.  ■  c  .  T  =cT  and  the  sum 

K  pK  X  p  X 

I  hk(pyk)x  in  Eq.  (3)  vanishes,  since  j  YkVk  *  0. 


The  boundary  conditions  of  the  unbumt  mixture'  have  been  specified  as: 

T  =«  298. 15K  for  a  stoichiometric  hydrogen-air  mixture.  Since  the  overall  chemi¬ 
cal  equation  can  be  written 


2H2  ♦  02  +  N2  =  2H20  +  N2 

and  since  air  is  composed5  of  79%  N2  and  21%  02,  we  find  the  mole  (mass) 
fractions  at  the  unburnt  boundary  to  be:  Xq  =  0.1479  (0.2264),  X^  =  0.29S8 
(0.0283)  and  X  *  0.5563  (0.7453).  2 

N2 

The  boundary  conditions  of  the  burnt  mixture  are 


Tx  •  (Yk)x  -  0  (k=l,2,...)»  (12) 

The  Adiabatic  flame  temp,  Tg,  is  2383°K  as  computed  by  the  NASA-Lewis  Codeb. 

To  define  the  starting  profiles  we  let  ip.  and  be  the  left  and  right 
boundaries  of  the  integral  of  integration.  Then  define? 

-  <Pt  +  0.24  (*R  -  <PL)-  (13a) 

^2  *  +l  +  0,64  *  (13b) 


The  composition  of  dry  air  is  by  mole  X;  78.084  R2*  20.946  0%,  0.934  Ar, 

0.033  CO2,  0.003  rare  gases.  See  Matheson  Gas  Data  Book ,  Sth  Ed.  by  W.  Broker 
and  A.L.  Mossman ,  Matheson  Gas  Products ,  E.  Rutherford,  N.J.  1971 ,  p.  9. 

ft 

Gordon,  S.  and  McBride,  B.J.,  Computer  Program  for  Calculation  of  Complex 
Chemical  Equilibrium  Compositions ,  Rocket  Performance,  Incident  and  Reflected 
Shocks  and  Chapman-Jouquet  Detonations,  NASA-SP-273,  1971,  (1976  program 
version) . 

7 Coffee,  T.P.  and  Heimerl,  J.M.,  BRL  Technical  Report  ARBRL-TR-0221 2,  Jan.  1980. 
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and 


v«  ■ 


Yku  •  *li'i  *1 
Yk„  *  -  Yku>  sl"  <7  < 

Yk,  •  *2S  *i*R- 


.  .  *-  *1  ,2,2  ,  ,  , 
ipjf3  >  -  *1  i  *  i  *2 


(13c) 


T  is  defined  similarly.  The  choice  of  the  particular  function  that  defines 
the  Yk  and  T  between  ipi  and  ij>2  is  not  important.  We  have  used  a  straight  line 
with  success,  but  defining  a  smooth  function  is  slightly  more  efficient. 

III.  The  Numerical  Solution 

g 

A  modification  of  the  PDECOL  package  was  used  to  solve  this  problem  . 
PDECOL,  a  rather  general  package  for  solving  one-dimensional  differential  equ¬ 
ations,  employs  the  method  of  lines  in  finding  a  solution. 


Spatial  discretization  is  accomplished  by  finite  element  collocation 
methods  based  on  B-splines.  The  basic  assumption*  is  that  the  solution  can  be 
written  in  the  form 


and 


NC  ... 

Y  =  I  CUj  (t)  B.  (i/i) ,  k  =  1 

k  i=i  K  i 

_  NC  (i) 

T  =  Z  C  (t)  B.  (*) , 
i=l  N+l 


.,N 


(14a) 

(14b) 


where  the  basis  functions  Bj^  (iji),  i=l,...,NC  are  B-splines  and  span  the  solu¬ 
tion  space  for  any  fixed  t.  to  within  a  small  error  tolerance.  The  time- 
dependent  coefficients  C^(i)  are  determined  uniquely  by  requiring  that  the 
preceding  expansion  satisfy  Eqs.  (6)  and  (7)  exactly  at  (NC  -2)  interior 
(collocation)  points.  The  boundary  conditions,  Eq.  (4),  are  also  converted  in¬ 
to  ordinary  differential  equations.  If  there  is  a  null  boundary  condition,  an 
extra  collocation  point  is  added.  The  &i  are  piecewise  polynomials  of  order 
KORD.  The  user  supplies  a  set  of  NB  breakpoints  (i.e.  a  set  of  strictly  in¬ 
creasing  locations  where  the  polynomials  are  joined)  and  the  value  for  NCC, 
the  number  of  continuity  conditions  to  be  applied  at  the  breakpoints.  From 
this  information  PDECOL  generates  a  set:  NC  =  KORD(NB-l)  -  NCC(NB-2)  basis 
functions  and  collocation  points.  Since  by  definition  a  B-spline  is  zero 
except  over  a  small  interval,  at  any  collocation  point  no  more  than  KORD  of 
the  B-splines  are  nonzero.  So  the  system  of  ODE’s  for  the  coefficients  (i3 
is  not  fully  coupled. 

This  system  is  integrated  in  time  by  using  a  variant  of  the  Gear  integrator 
It  is  at  this  point  in  the  solution  that  we  have  modified  the  PDECOL  package**. 
Normally  all  the  Yfc's  and  T  would  be  solved  simultaneously.  This  procedure 
would  consume  inordinate  amounts  of  computer  storage.  Since  we  are  interested 

TJ - 

Coffee ,  T.P.,  BRL  Memorandum  Report  to  be  published.  ' 


r 


12 


only  in  steady-state  solutions,  we  have  adapted  the  method  of  successive  cal¬ 
culation.  At  each  time  step  this  process  considers  each  PDE  in  turn  and  values 
for  other  dependent  variables  are  taken  from  the  most  recent  computation.  In 
the  steady-state  this  procedure  is  justified  because  0^)^.  =  Tt  =  0  and  the 
system  is  in  fact  decoupled  in  time. 

The  code  requires  the  selection  of  the  following  parameters.  The  thermal, 
ip-space  and  temporal  normalization  factors  were  taken  to  be  1000K,  5x10“^  gm  cm- 2 
and  lxl0-3  sec,  respectively.  ij>0,  the  cold  boundary  location  is  taken  to  be 
zero  while  1J/5  the  normalized  "boundary  at  infinity"  is  taken  to  be  2.90. 

The  time  integration  is  controlled  by  a  user  supplied  error  tolerance  e. 
Single  step  error  estimates  divided  by  CMAXfcCi)  will  be  kept  less  than  c  in 
the  root-mean-square  norm.  In  the  present  case  e  was  chosen  as  3xl0-^.  In 
PDECOL,  CMAXfcCi)  is  initially  set  to  the  maximum  of  | c^ Ci Y |  and  1.0.  There¬ 
after,  CMAXfc (i)  is  the  largest  value  of  |cfc(i)|seen  so  far,  or  the  initial 
CMAX^(i)  if  that  is  larger.  This  error  criterion  is  not  appropriate  for  flames. 

In  a  flame,  radical  species  with  relatively  small  concentrations  will 
control  the  flame,  and  we  want  to  compute  these  accurately.  However,  some 
species  concentrations  at  some  locations  will  approach  zero,  and  we  do  not 
want  to  waste  time  computing  a  negligible  concentration  very  accurately. 

So  we  will  use  a  semi-relative  error  control.  CMAX^fi)  will  be  chosen 
as  the  maximum  of  c]<(i)  and  a  user  supplied  parameter  SREC.  So  concentrations 
less  than  SREC  will  be  computed  less  accurately.  Here  we  have  used  SREC  =  10-6. 

The  major  difficulty  in  efficiently  solving  the  flame  equations  is  choosing 
an  appropriate  set  of  breakpoints.  These  must  be  close  enough  that  spatial 
errors  do  not  destroy  the  solution  yet  not  so  dense  that  one's  computer  resources 
are  exceeded.  The  breakpoints  should  be  densest  in  the  flame  front,  where  the 
gradients  are  very  steep. 

This  problem  can  be  dealt  with  by  working  with  mQ,  the  mass  flux  through 
the  origin.  This  can  be  iteratively  modified  to  match  the  mass  flux  through 
the  flame.  It  can  also  be  modified  to  move  the  flame  front  towards  the  center 
of  the  interval  of  integration,  if  necessary.  The  details  of  this  procedure 
are  discussed  in  Reference  7. 

So  we  may  choose  our  breakpoints  to  be  densest  in  the  center  of  our 
interval  of  integration,  we  have  developed  a  procedure  to  generate  this  type 
of  breakpoint  sequence. 

The  user  must  supply  NINT,  the  number  of  intervals  CNB  =  NINT  '  1),  NCN, 
the  number  of  intervals  of  equal  length  that  will  be  the  center  of  the 
interval,  and  FC,  the  ratio  between  the  longest  intervals  (on  the  boundaries) 
and  the  shortest  intervals.  Also  let  L  be  the  total  length  of  the  interval  of 
integration. 

The  program  generates  a  set  of  intervals  whose  lengths  increase  by  a 
constant  factor  a,  where 

a  =  log-1  (2(log  FC)/ (NINT  -  NCN)).  (15a) 
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The  common  length  LC  of  the  NCN  center  shortest  intervals  is 

LC  =  L/ (NCN  +  2a  (a(NINT"NCNi/2-l)/Ca-l)) .  ( 

Some  experimentation  is  necessary  to  choose  the  proper  values  of  the 
above  parameters.  In  the  present  case  we  have  chosen  NINT  =  24,  NCN  =  6  and 
FC  =  6.  Finally  the  order  of  the  B-splines  KORD  is  taken  to  be  4  and  the 
number  of  continuity  conditions  at  the  collocation  points,  NCC,  is  set  equal 
to  2. 


The  code  was  executed  on  the  BRL  CYBER  76  system:  it  required  approxi¬ 
mately  1220008  words  small  core  memory  and  took  about  100  seconds  to  run  from 
a  restart  file. 


IV.  Input  Coefficients 


Kinetics 


The  steady-state,  stoichiometric,  plane  hydrogen-air  flame  defined  by 
Peters  has  the  kinetic  parameters  given  in  Table  1.  (k  =  BT11  exp(c/T)  in  units 
of  cm,  mol,  s,  K.)  As  a  check  on  the  operation  of  the  code  we  employed 
another  kinetic  network  listed  by  Dixon-Lewis^.  These  parameters  are  given  in 
Table  2.  We  remind  the  reader  that,  at  this  time.  Table  2  contains  our  best 
understanding  of  the  H2/O2/N2  system  and  that  Table  1  is  a  simplified  set  for 
the  purposes  of  this  exercise. 

Thermodynamics 

The  enthalpy  and  heat  capacity  are  functions  of  species  and  of  temperature. 
The  JANAF  values  for  the  eight  species  H,  OH,  0,  HO2,  H2,  O2,  H2O  and  N2  are 
adequately  described  by  the  polynomial  fits  of  Gordon  and  McBride6.  These  fits 
take  the  form 

5 

Ej.0  =  R  (aQ  +  E  an  T^/n)  (16a) 

n=l 

and  g 

C  =  R  E  a  Tn_1  (16b) 

p  .  n  J 

c  n=l 

where  the  coefficients  are  listed  in  Table  3  The  heat  capacity  and  specific 
heat  at  2000K  are  listed  in  Table  4  and  were  requested  for  use  at  the  workshop. 

When  the  values  from  these  polynomials  are  compared  to  the  tabular 
valuesl°-12>  those  for  HO2  agree  to  better  than  3.6%  while  all  the  others  agree 
to  better  than  1%. 

0 

Dixon-Lewis,  G. ,  Philos.  Trans.  Roy.  Soo.  (London)  292,  45-99,  1979. 

10JAMF  Thermo chemiaal  Tables ,  2nd  Ed.  by  D.R.  Stull  and  H.  Prophet,  NSRDS-NB37 , 
June  1971. 

11 

Chase,  M.W.,  Curnutt,  J.L.,  Hu  ,  A.J.,  Prophet,  H .,  Syverud,  A.N.,  and  Walker, 
L.C.,  J.  Chem.  Ref.  Data  3,  311-480,  1974. 

12 

Chase,  M.W.,  Cumutt,  J.L.,  Prophet,  H.,  McDonald ,  R.A.,  and  Syverud,  A.N.,  J. 
Phys.  Chem.  Ref.  Data,  4,  1-175,1975. 
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Table  1.  Peters*  l^-Air  Kinetics 


# 

Reaction 

B 

n 

c 

Ref. 

1 

H2  +  OH  =  H20  +  H 

2.2(13) 

0 

-2590 

A.p77 

2 

02  +  H  =  0  +  OH 

2.2(14) 

0 

-L  ,50 

A,p9 

3 

H2  +  0  =  H  +  OH 

1.8(10) 

1.0 

-4480 

A,p49 

4 

0H+0H  =  H20  +  0 

6.3(12) 

0 

-550 

A,pll9 

5. 

H+H+M  =  H2+M 

6.4(17) 

-1 

0 

A,p299* 

6. 

OH  +  0+M  =  H02+M 

5.  (16)*** 

0 

0 

A,p411 

7. 

o2+h  +  m  =  ho2  +  m 

1.5(15) 

0 

+550 

A,p377* 

8. 

H02+H  =  OH  +  OH 

2.5(14) 

0 

-950 

A,pl45 

9. 

H02  +  H  =  H2+02 

2.5(13) 

0 

-350 

A,pl61 

10. 

OH  H02  *  H20+02 

1.5(13)**** 

0 

0 

A,p251 

11. 

H02  +  0  =  02  +  OH 

6.3  (13) 

0 

-350 

B,pl73 

Values  listed  are  for  argon  as 

third  body. 

**  Values  listed  are  for  helium  or  argon  as  tmrd  body. 
***  Value  of  Petersen  $  Kretschmer  at  300#K  ;  for  M  =  02 
****  Appears  estimated  based  on  Leeds  values. 


A  -  Evaluated  Kinetic  Data  For  High  Temperature  Reactions  by  D.L.  Baulch, 

D.D.  Drysdale,  D.G.  Home,  A.C.L.  Lloyd,  Vol.  I,  CRC  Press,  Cleveland,  1972 

B  -  Evaluated  Kinetic  Data  for  High  Temperature  Reactions,  by  D.L.  Baulch, 

D.D.  Drysdale,  J.  Duxbury  and  S.J.  Grant,  Vol.  3,  Butterworth,  Boston,  1976 
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Table  2.  The  Kinetic  Scheme  of  Dixon- Lewis 


# 

Reaction 

B 

n 

C_ 

1 

OH  +  H2  =  H20+H 

1.17(9) 

1.3 

-1825 

2 

H  +  02  =  0H+0 

1.42(14) 

0.0 

-8250 

3 

0  +  H2  =  OH+H 

1.8  (10) 

1.0 

-4480 

4 

h+o2+h2  *  ho2+h2c 

1.03(18) 

-0.72 

0.0 

5 

H+H02  *  OH+OH 

1.4(14) 

0.0 

-  540 

6 

h+ho2  *  o+h2o 

1.0(13) 

0.0 

-  540 

7 

H+H02  *  H2+02 

1.25(13) 

0.0 

0.0 

8 

0H+H02  =  H20+02 

7.5(12) 

0.0 

0.0 

9 

0+H02  a  OH+02a 

1.4(13) 

0.0 

-  540 

10 

0+H02  a  0H+02a 

1.25(12) 

0.0 

0.0 

11 

H+H+H2  =  H2+H2 

9.2(16) 

-0.6 

0.0 

12 

H+H+N2  a  H2+N2 

1.0(18) 

-1.0 

•  0.0 

13 

H+H+Oz  =  H2+02 

1.0(18) 

-1.0 

0.0 

14 

H+H+H20  =  H2+H20 

6.0(19) 

-1.25 

0.0 

15 

H+QH+M  a  HjO+M^ 

1.6(22) 

-2.0 

0.0 

16 

H+O+M  »  OH+^ 

6.2(16) 

-0.6 

0.0 

17 

OH+OH  =  0  ♦  H20 

5.75(12) 

0.0 

-  390 

a, 

Dixon-Lewis  defines  the  rate  coefficient  for  0  ♦  H02  *  OH  +  as  +  k1Q» 
bValue  for  M*H2,N2  or  02;  X5  for  M=H20 . 

cChaperon  efficiencies  relative  to  H_  =  1.0  are  0.44,  0.35  and  6.5 
for  N2,  02  and  HjO,  respectively.  1 


16 


TABLE  3.  THERMODYNAMIC  COEFFICIENTS 


Table  4.  Heat  Capacity  and  Specific  Heat  at  2000K. 

C_(cal/mole-K)  c  (Joule/ ga-K) 

Species  — - -  -£ - - - 


H 

4.968 

20.79 

OH 

8.275 

1.923 

0 

4.976 

1.301 

ho2 

12.906 

1.636 

H2 

8.173 

17.10 

°2 

9.045 

1.183 

h2o 

12.199 

2.8356 

N2 

8.609 

1.286 
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TRANSPORT 


Individual  Species 

There  are  two  classes  of  transport  parameters; Xi,  the  heat  capacity  of 
the  ith  species,  and  Djj,  the  binary  diffusion  coefficient.  We  consider  them 
in  turn . 


The  theoretical  expression  describing  the  thermal  conductivity  of  hard 
sphere  (i.e.,  atomic  species)  is  simply  related  to  the  viscosity;  specifically. 


Wni 


15R/4. 


To  account  for  the  fact  that  energy  storage  in  a  polyatomic  molecule  can  take 
place  through  inodes  other  than  translational  we  have  examined  the  Eucken  and 
modified  Eucken  correlations.  In  general  form, 

X.M./n.  *  f„  CL  +  f.  S-  *  •  (18) 

1  1  i  tr  tr  xnt  int,i  v  1 

where  ftr  is  set  equal  to  5/2  so  that  this  expression  reduces  to  the  monatomic 
case  (i.e.  when  f^nt  =  0,  there  Ctr  =  3/2R) .  Cint  ^  is  taken  to  be  Cv  i-Ctr  and 
there  remains  the  evaluation  of  fint.  Eucken  set  fjnt  to  unity.  Others  have 
suggested  that  the  flow  of  internal  energy  is  a  type  of  diffusional  transport 
and  have  found  that  fint  »  (Sc)*l,  where  the  Schmidt  number.  Sc,  is  n/pDs.  Ds 
the  self  diffusion  coefficient,  *  1.20  pnJl(^»2)/n(l ,1)  (here  the  ideal  gas 
has  been  assumed).  Since  the  ratio  of  collision  integrals  is  approximately 
1.1  and  is  almost  independent  of  temperature,  we  find  (Sc)"l  *  1.32.  This 
value  (f^t  *  1.32)  leads  us  to  the  modified  Eucken  correlation. 

Usually  experimental  values  of  X^  for  non-polar  gases  lie  between  the 
values  calculated  by  the  two  Eucken  forms  and  so  we  have  elected  to  take  the 
arithmetic  mean  of  the  two  forms. 

The  ensuing  expression  reduces  to 


x 


1.20R  +  2.32  Cpi 


where  we  have  used  the  relation  R  =  Cp-Cv.  C  can  be  computed  from  the 
thermodynamic  coefficients  of  Gordon  and  McBride^,  Eq.  16b. 

We  see  from  Eq.  19  that  to  find  the  expression  for  Xi»  we  need  to  find  the 
expression  for  ni- 

The  viscosity  of  the  ith  species  is  given  by: 

r,i(gm-cm"1s'1)  =  2.699  x  10"5  (K^T)  /ojV2,25  (t/)  (20) 

The  discussion  follows  that  of  Reid,  R.C.  and  Sherwood,  T.K.,  The  Properties 
of  Gases  and  Liquids  ,  2nd  Ed.  MoGraw-Hill ,  N.T.,  1966.  Viscosity  Chap.  9, 
Thermal  Conductivity ,  Chap.  10  and  Diffusion  Coefficients ,  Chap.  11. 
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*  -1 

where  Tj  s  (e/k) i  T.  The  Lennard- Jones  (or  Stockmayer,  in  the  case  of  H2O) 
parameters  are  given  in  Table  5. 


TABLE  5.  LENNARD- JONES/STOCKMAYER  PARAMETERS 


Species 

a 

(10'8cm) 

TO 

Remarks 

H 

2.05 

145 

a 

OH 

2.947 

127.2 

like  0 

0 

2.947 

127.2 

b 

ho2 

3.372 

128.7 

like  02 

H2 

2.92 

38.0 

a 

°2 

3.372 

128.7 

b 

h2° 

2.60 

572 

a 

3.62 

97.5 

a 

aJ.  Wamatz  "Calculation  of  the  Structure  of  Laminar  Flat  Flames  II:  Flame 
Velocity  and  Structure  of  Freely  Propagating  Hydrogen-Oxygen  and  Hydrogen 
Air-Flames"  Ber.  Bunsengs  Phys.  Chem.  82,  643-649  (1978). 

kj.M.  Heimerl  and  T.P.  Coffee,  "The  Detailed  Modeling  of  Premixed  Laminar 
Steady-State  Flames.  I  Ozone."  Combustion  and  Flame  39,  301-315  (1980). 


In  order  to  evaluate  the  collision  integrals  for  polar  molecules  not  only 
is  e/k  required  but  also  the  dimensionless  parameter 

6  *  mH20^  eH20  °H2oj'  (21) 

The  dipole  moment  for  H2O,  vh20»  is  taken  to  be  1.844  Debye  (l  Debye  =  101® 
dyne4-cm2).  This  value  agrees  with  those  reported  by  McCellanl*  (1-82)  and 
Nelson*5  et  al.  (1.85).  With  this  value  for  and  the  parameters  from  Table 
4  we  find  6  =  1.225. 


14 

MaCellan ,  A.L. ,  Tables  of  Experimental  Dipole  Moments ,  Freeman,  San  Franoieoo, 
1963. 

*^Nel8on,  R.D.,  Jr.3  hide,  D.R.,  Jr.  and  Maryott,  A. A.,  NSRDS-iiBSlO 3  Sept.  1967. 
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We  have  found  it  computationally  efficient  to  use  an  analytic  fit  for 
X.  and  ni-  Values  for  these  parameters  were  generated  from  equations  (19)  and 
(20),  respectively.  The  temperature  range  was  300-3000K  and  evaluations  were 
made  at  100K  intervals.  The  functional  form  AT®  was  then  fitted  to  these  data 
in  a  least  squares  sense.  The  coefficients  A  and  B  for  each  species  and  the 
overall  rms  error  of  each  fit  is  listed  in  Table  6.  Table  7  lists  the  species 
thermal  conductivities  at  2000K. 

13 

The  binary  diffusion  coefficients  are  generated  from  the  expression 
pD„  (cm2s-1atmos-1)  =  1.858  x  10-3T3/2 [ (Mj+M^. )  (M^)"1]*/^?^1,15  (T*.) 

(22) 

where  T?.  =  (e/k). .  T.  For  non-polar,  non-polar  interactions  the  combining 
laws  are  e„  =  (e-e..)1  (23) 

and 

aij  =  i(oi  +°j)*  (24) 


For  polar,  non-polar  interactions  we  employ  the  modified  combining  laws: 
a 


np 


and 


'np 


<«„  *  v  5-1/6 

(25) 

(26) 

♦ 2-1  (vy.!) 

lep/En]i 

(27) 

8"^  M  2/e  o  3. 

P  P  P 

(28) 

where 


and 

t  *  = 

P 

We  note  in  passing  that  26  =  t  *.  To  evaluate  equation  (27)  we  need  to  specify 
the  polarizability  of  the  non-polar  molecule,  a  .  For  H2,  02  and  N2,  a_  (in  units 
of  10-24  cm3)  iglo.  0.79,  1.60  and  1.76,  respSctively. 

The  values  of  the  Djj  are  computed  at  100K  intervals  over  the  temperature 
range  300-3000K  by  means  of  Eq.  (22)  and,  as  the  values  for  Xi,  are  fitted  to 
the  convenient  functional  form  AT&.  The  coefficients  A  and  B  and  the  overall 
error  of  the  fit  (in  an  rms  sense)  are  given  in  Table  8. 

Some  of  these  values  can  be  checked  against  experimental  data  and  other 
Dij  determinations.  Table  9  shows  values  of  Dij  computed  from  the  parameters 
listed  in  Table  8  and  experimental  data  taken  at  273K  and  one  atmosphere  pres¬ 
sure.  We  find  this  agreement  good  to  excel lant. 


TABLE  7.  THERMAL  COWUCTIVITY  AT  2000X. 
(computed  from  Table  6) 


Xi 

Xi 

Species 

cal-cm" 

-K-1s'1 

(W/mKj 

H 

2.6822 

* 

>-* 

O 

1 

Ol 

1.122 

OH 

4.847 

* 

1 

O 

H 

0.2028 

0 

3.318 

* 

10-4 

0.1388 

ho2 

3.900 

* 

10-4 

0.1632 

H2 

1.712 

* 

to 

1 

O 

H 

0.7163 

°2 

2.921 

* 

O 

1 

0.1222 

H2° 

5.892 

* 

10-4 

0.2465 

N, 

2.687 

* 

* 

O 

1 

0.1124 

Table  10  shows  the  expressions  of  Marrero  and  Mason1  for  selected  D. . , 
also  in  the  form  AT8.  Their  expressions  are  based  upon  data  (except  for  J 
H2-N2)  and  they  rate  the  quality  of  their  fits  in  the  following  order: 


h2°  -  n2 

±  4% 

H2°  “  °2 

±  7% 

H  -  H2 

i  5%  (300°K) ;  ±  30%  (>1000K) 

h-n2,  o-n 

2,  0-02  ±  10% (300K) ;  ±  25%(>1000K). 

The  last  column  of  Table  10  shows  the  comparison  of  their  values  with  the 
ones  generated  from  the  parameters  listed  in  Table  8.  This  comparison  is  made 
over  the  interval  300  £T<_ 3000K  except  where  noted.  Comparisons  for  the  major 
species  are  within  13%  (rms) .  The  O-N2  and  O-O2  show  larger  rms  deviations 
but  their  effects  on  flame  speed  and  profiles  is  expected  to  be  small. 

- - 

Marrero ,  T.R.  and  Mason,  E.A.,  J.Phja.  Chem.  Ref .  Data  2,  3-118,  1972. 
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Table  8.  ANALYTIC  FIT  OF  FORM  AT8  TO  THE  BINARY  DIFFUSION  COEFFICIENTS  D„ 


pD.  . 

A 

B 

Error 

i 

i 

(X 1 0“ ^ cm^  s  ”  *  atm "" 1 ) 

(%) 

H 

H 

2.1700 

1.6846 

1.3 

H 

OH 

1.0997 

1.6815 

1.2 

H 

0 

1.1016 

1.6815 

1.2 

H 

ho2 

0.91865 

1.6818 

1.2 

H 

H2 

1.6560 

1.6649 

0.4 

H 

°2 

0.91907 

1.6818 

1.2 

H 

H2° 

0.6409 

1.7525 

2.7 

H 

N2 

0.90158 

1.6758 

0.9 

OH 

oh 

0.27234 

1.6786 

1.1 

OH 

0 

0.27656 

1.6786 

1.1 

OH 

ho2 

0.20564 

1.6788 

1.1 

OH 

»2 

0.73760 

1.6640 

0.3 

OH 

°2 

0.20673 

1.6788 

1.1 

OH 

h2o 

0.16361 

1.7437 

2.6 

OH 

N2 

0.20894 

1.6734 

0.8 

0 

0 

0.28072 

1.6786 

1.1 

0 

H°2 

0.20984 

1.6788 

1.1 

0 

»2 

0.74002 

1.664 

0.3 

0 

°2 

0.21091 

1.6788 

1.1 

0 

h2o 

0.16622 

1.7437 

2.6 

0 

N2 

0.21296 

1.6734 

0.8 

ho2 

ho2 

0.14846 

1.6791 

1.1 

ho2 

H2 

0.62386 

1.6641 

0.3 

H02 

°2 

0.14961 

1.6791 

1.1 

ho2 

h2o 

0.12144 

1.7445 

2.6 
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Table  8.  ANALYTIC  FIT  OF  FORM  AT6  TO  THE  BINARY  DIFFUSION  COEFFICIENTS  D.  ^  (cont'd) 


P°ii 

A 

B 

Error 

i 

_ _ 

i 

(X10"4cm2s'1atm'1) 

(%) 

ho2 

N2 

0.15364 

1.6736 

0.8 

H2 

H2 

1.1299 

1.6600 

0.2 

H2 

°2 

0.62442 

1.6641 

0.3 

H2 

h2o 

0.58694 

1.6937 

1.7 

H2 

N2 

0.59990 

1.6627 

0.3 

°2 

°2 

0.15076 

1.6791 

1.1 

°2 

h2o 

0.10965 

1.758S 

2.8 

°2 

N2 

0.15474 

1.6736 

0.8 

H2° 

h2o 

0.026558 

1.9523 

3.1 

h2o 

N2 

0.12198 

1.7403 

2.5 

N2 

N2 

0.15678 

1.6700 

0.6 

Note:  The  logarithmic  weighting  of  these  fits  was  taken  into  account.  See 
h.J.  Cvetanovic  and  D.L.  Singleton,  Int.,  J.  Chem.  Kinet.  9,  481-438, 
1977  anrl  9,  1007-1009,  1977. 
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TABLE  9.  COMPARISON  OF  THE  VALUES  OF  SELECTED  D...  FROM  TABLE  8  COMPUTED  AT 
1  atmos.  and  273K  WITH  EXPERIMENTAL  DATA 


Species  pair 

atb 

D  * 

exp 

H2 

-  N2 

0.674 

0.674 

H2 

-  h2o 

0.784 

0.759 

°2 

-  h20 

0.211 

0.214 

H2 

-  °2 

0.750 

0.701 

N2 

-  H20 

0.212 

0.212 

N2 

-  °2 

0.185 

0.181 

H  - 

H2 

1.88 

1.86 

0  - 

°2 

0.259 

0.281 

H  - 

N2 

1.09 

0.945 

0  - 

N2 

0.254 

0281 

H2 

-»2 

1.25 

1.26 

°2 

-  °2 

0.186 

0.184 

N2 

‘  N2 

0.184 

0.180 

*Dexp  =  273 *  P  =  1  atmos)  from:  Warnatz,  J. 

Ber  Bunsenges.  Phys  Chem  82  643-649,  1978,  Table  2. 
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TABLE  10.  COMPARISON  OF  SELECTED  D.^  WITH  THOSE  OF  MARRERO  AND  MASON17 


D.  . 

5  2 

10  A(atmo-cm  ) 

B 

Temperature 
Range  (K) 

Comparison  With 

Values  from  Table  8 

i 

i 

H2 

°2 

4.17 

1.732 

252-104 

9.8% 

N2 

°2 

1.13 

1.724 

285-104 

6.0% 

N2 

h2o 

0.187 

2.072 

282-373 

1.7%  (300K  only) 

{  0.189 

2.072 

282-450 

12.9%  ( 300-500 j 

°2 

h20 

l  2.78 

1.632 

450-1070 

10.6%  (500-1100) 

H 

H2 

11.3 

1.728 

190-104 

8.5% 

0 

N2 

1.32 

1.774 

280-104 

25.7% 

0 

°2 

1.32 

1.774 

280-104 

22.8% 

H„ 

No 

* 

* 

4 

65-10 

12.2% 

*«,n(pDH  N  )  =  in  A  +  s  £nT  -  £n[£n(if>o/kT)]2-(S/T)  -  (S'/T2) 
where 

103A  =  15.39  atraos-cin2  [s(K)S]  1 

s  =  1.548 

10'8<J>  /k  =  0.316 
o 

S  =  -  2.80 
S’  =  1067. 
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The  LJ  parameters  for  Dh2-N2  were  determined  from  the  empirical  combining 
rules  Eqs.  (23)  and  (24)  for  h2  and  N.,.17 

Table  11  shows  the  results  of  comparisons  of  the  viscosity  and  thermal 
conductivity  as  computed  from  the  parameters  on  Table  6  and  the  data  listed 
by  Touloukian  et  al'S,19.  The  large  deviation  for  A^  ^  suggests  that  a  fit 
could  be  made  to  the  data  themselves.  n2J 


TABLE  11. 

COMPARISON  OF  THERMAL  CONDUCTIVITY 
Touloukian  et  al. 18,19 

AND  VISCOSITY 

WITH  DATA  OF 

Species 

H2 

A  deviations 
+  4 . 7i 

Temperature 
Range  (K) 
300-2000 

n  deviations 
4.90 

Temperature 
Range  (K) 
300-2000 

N2 

-  3.6% 

300-3000 

1.34 

300-2000 

°2 

+  1.1% 

300-1500 

3.03 

300-2000 

h2o 

-  29.8% 

300-900 

4.13 

300-1000 

+  =  data  >  computation 
Multispecies  Transport 

Once  the  and  D^-  are  determined,  as  above,  there  remains  to  assemble 
the  thermal  conductivity  of  the  mixture  and  decide  how  to  handle  the  multi¬ 
species  diffusion. 

18 

Although  some  rather  involved  formulations  are  available  in  the  literature 
to  determine  Amix, we  have  found2^  that  the  values  from  these  expressions  differ 
little  from  the  expression 


A  . 


mix 


0.5  (E.  X.  A .  +  (E  X./A.)"1}  . 
i  1.  i  vi  i 


(29) 


For  the  diffusion  transport  parameters  we  employ  a  generalization  of 
Fick's  Law,  namely 

JO 

Touloukian ,  7.5.,  Liley,  P.E.  and  Saxena ,  5.(7.,  Thermal  Conductivity  (non- 
metallia  Liquids  and  Gases),  IFI/Plenvm,  N. 7. ,  1970. 

19 

Touloukian ,  7.5.,  Saxena ,  S.C.  and  Hestermans ,  P.  Viscosity .  IFI/Plertum ,  N.  7. , 
1975. 


pi 9 

Coffee ,  T.P.  and  Heimerl,  J.M.t  Combustion  and  Flame  43j  273-289  (1981)  and  BRL 
Tech.  Report  ARBRL-TR~02302t  Mar.  81. 
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This  is  the  so-called  Method  V  considered  in  Ref.  20.  Method  VI  is  considered 
in  the  Appendix. 
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(30) 


Vi 


"  Dim  <Vx 


where 


D. 

lm 


(1-X. ) (I .  , . X./D.  .) 
1  jfi  r  i3J 


-1 


(31) 
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Again  we  have  found  little  difference  in  the  computed  flame  speeds,  species 
profiles  or  temperature  profiles  using  this  generalized  Fick's  law  vs.  more 
complex  formulations  that  solve  for  the  diffusion  velocity  (such  as  the  Stefan- 
Maxwell  equations  do) . 


We  have  further  noted  that  in  the  Lennard-Jones  formalism 


D.  .  «  T1*5/^1,1)  , 
ij 


(32) 


Now  n^1,1^  «  T  and  since  p  =  T  *  we  see  that  p2D„  is  nearly 

independent  of  temperature.  Generalizing  we  assume  that  P^D|m  is  constant. 

Likewise  for  a  monatomic  gas  A«  T*/V2,V  where  $2^2*2^  «  -0.16.  It  is  thus 

reasonable  to  assume  that  pX  is  independent  of  temperature. 

We  now  outline  a  procedure  that  permits  an  a  priori  selection  of  these 
quantities. 


For  a  given  flame  we  know  Tu  and  Yiu>  the  temperature  and  the  mass  fraction 
of  the  unburnt  mixture.  We  also  know  the  kinetic  scheme,  the  specific  heats 
and  the  specific  enthalpies.  Since  enthalpy  is  conserved,  i.e.. 


N 

I 

i=l 


Y.h. 

i  i. 


N 

l 

i=l 


Y.h. 

l  l 


(33) 


we  can  compute  the  adiabatic  flame  temperature  Tg,  as  well  as  the  mass  fractions 
of  the  burnt  mixture,  Yig.  The  numerical  procedure  is  to  guess  a 
trial  temperature  and,  using  an  ODE  package,  find  the  corresponding  equili¬ 
brium  mass  fractions  Yjt*  We  then  compute  the  enthalpy  of  the  burnt  mixture 
and  compare  with  the  unburnt  mixture  enthalpy.  Tt  is  iteratively  adjusted 
until  the  burned  and  unburned  mixture  enthalpies  agree  to  within  a  predeter¬ 
mined  error  tolerance.  We  then  accept  these  values  of  T.  and  Y.„ 

As  a  heuristic  rule  we  select 

T  =  0.5(TB  -  Tu)  (34a) 


and 


Y. 

l 


0.5  (Y.b  ♦ 


Y.  ) 
iu' 


(34b) 
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then  evaluate  Dim  and  \i  from  Eqs.  (31)  and  (29)  respectively.  Then  p  Dim, 
i  *  1.2...N-1  andpXi  are  evaluated.  The  diffusion  velocity  is  found  from 
the  constraint  on  the  diffusive  mass  flux  that 


lY^i 


(35) 


At  this  level  of  approximation  we  find  it  convenient  to  employ  a  constant 
or  global  Cp.  So  we  assume  cn.>  cn  *  constant  for  all  i.  Then  the  relation 


».  =  h  °  *  f 
1  1  Jl 


pi  p 


c  .  dT 
T  pl 

*  rv 


(36) 


becomes 


h.  -  h.  ♦  c  (T  -  T  )  and 
11  p  o 


(37) 


the  mixture  enthalpy  is  given  by 


l  Y.h.  =  E  Y.h.  ♦  c  (T-T  ) 
11  ii  p  o 


(38) 


Substituting  this  expression  into  Eq.  (33)  we  find 


CP 


*  hi°tYiu 


YiB> 


Td-T 
B  u 


(39) 


The  values  used  here  are  found  in  Table  12  and  yield  an  average  Cp  = 
3.65  x  10-1  cal-mole-lK-1.  we  also  find  that  p(T  =  1344K)  =  2.0364  x  10'4 

2390K.  This  value  compares  well  with  the  HASA-Lewis 


gm  cmg3 


and  that  T, 


value  of  2383K. 


B 


Table  13  shows  the  Di,  evaluated  at  T  =  2000K  as  required  for  use  in  the 
workshop  comparisons;  the  values  at  the  mean  temperature,  T  =  1344K,  are  in¬ 
cluded  for  completeness. 


and  Dxm>  both  at  1344K,  used  in  the  compu- 


Table  14  contains  values  of  Xj  »uu  uim> 
tation  of  Dim  and  pXmix.  The  value  of  this  density  at  1344,  given  above, 
was  computed  from  p  *  p/(RT£Yi/Mi) ;  the  mass  fractions,  Yi,  are  taken  from 
Table  12  (ie.,  Y1344);  and  Xmix  is  composed  from  Eq.  (29).  This  equation  re¬ 
quires  the  species  mole  fractions  and  these  are  computed  from  the  standard 

relat  •  ¥i  =  fVj/Mii/r  £V •  /Mj  ^  Thft  ualnoe  f Ar 

code 
3.29  x 
is  6.70 


ion:  Xi  =  (Yi/Mi)/(^Yj/Mj) .  The  values  for  P*Dim  used  in  the  flame 
are  listed  in  the  last  column  of  Table  14.  The  value  found  for  Xmix  i 
x  10'4cal-s" *-cm  -IT1  and  |o  the  vp.lue  for  pX_.„  used  in  the  flame  co 


x  10 


t  2  1 

0  gm -cal -cm* -s'1 


mix 


is 

code 
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The  workshop  had  requested  values  for  at  2000K.  Table  15  contains 
a  listing  of  the  values  for  Yi  at  2000K  taken  from  the  flame  code  output.  Also 
listed  are  the  values  of  the  corresponding  mole  fractions.  This  information, 
the  values  of  Xi  from  Table  6  and  Eq.  (29)  were  used  to  complete  Xmix  * 

4.44  x  10-4  cal-s-l-K-1  fi.e.  0.1857  W  m-lK"1). 


TABLE  12.  PARAMETERS  NEEDED  TO  DETERMINE  DIFFUSION  TRANSPORT  PARAMETERS, 
EQUATIONS  (31)  and  (39) 


Y 

u 

y_b 

Y 

*1344 

^cal/gm) 

H 

* 

7 . 51 7 (-5) 

3.759(-5) 

5.21(+4) 

OH 

* 

5.081(-3) 

2.541 (-3) 

5. 58 (+2) 

0 

* 

4. 059 (-4) 

2. 030 (-4) 

3.72(+3) 

ho2 

* 

9.475(-7) 

4.738(-7) 

1.51 (+2) 

H2 

2 . 830(-2) 

1.212(-3) 

1.476(-2) 

-1.60(-1) 

°2 

2 . 264(-l) 

7.504(-3) 

1 . 170(-1) 

-3.57 (-2) 

h2o 

* 

2.404(-l) 

1.202(-1) 

-2.21(*3) 

N2 

7.4528C-1) 

7.453(-l) 

7.453(-l) 

-3.95(-2) 

T(°K) 

298 

2390 

1344. 

•Values  are  machine  zero. 


TABLE  13.  EVALUATION  OF 

D. .  AT  T  =  1344K  and 
i  J 

2000K 

i 

DL.(T  =  1344°K) 

D. . (T  =  2000 K) 

i 

*•3  _ 

\  -1 - 

cm  -s 

cm  -s 

H 

OH 

20.03 

39.08 

H 

0 

20.07 

39.15 

H 

HO 

16.77 

32.72 

H 

H 

26.76 

51.87 

H 

°2 

16.78 

32.74 

H 

H-0 

19.87 

39.87 

H 

n2 

15.76 

30.68 

OH 

0 

4.93 

9.61 

OH 

HO- 

3.67 

7.16 

OH 

h2 

11.84 

22.95 

OH 

°2 

3.69 

7.20 

OH 

h2o 

4.66 

9.33 

OH 

3.59 

6.98 

0 

h02 

3.74 

7.31 

0 

**2 

11.88 

23.02 

0 

°2 

37.68 

7.34 

0 

HO 

47.39 

9.48 

0 

«2 

36.59 

7.11 

ho2 

«2 

10.02 

19.42 

ho2 

°2 

2.68 

5.22 

HO- 

HO 

3.48 

6.97 

HO- 

"2 

2.64 

5.14 

°2 

10.03 

19.44 

h2 

H-0 

11.67 

22.88 

» 

n2 

9.54 

18.48 

02 

H-0 

3.48 

6.70 

°2 

»2 

2.66 

5.18 

H20 

N2 

3.39 

6.78 
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TABLE  14.  THE  TRANSPORT  PARAMETERS  AT  THE  MEAN  TEMPERATURES,  T-1344K 


Species 

Xi 

Di« 

p2D. 
y  ia 

X10”4cal-s"1-c»"1' 

-K"1 

xlO-'7  g»2-cni 

H 

20.52 

17.61 

7.30 

OH 

3.57 

4.24 

1.76 

0 

2.56 

4.32 

1.79 

ho2 

2.78 

3.15 

1.31 

H2 

12.67 

9.93 

4.12 

°2 

2.14 

3.21 

1.33 

h2o 

3.69 

3.95 

1.64 

N2 

1.98 

4.30 

1.78 

TABLE 

15.  FLAME  CODE  VALUES  OF  Y.^ 

AND  X.  AT  T=2000K 

Species 

h 

H 

1 

7. 79 (-4) 

1.84(-2) 

OH 

17 

1.00(-2) 

1.39 (-2) 

0 

16 

3. 54 (-3) 

5.22(-3) 

ho2 

33 

9.37(-6) 

6.70(-6) 

H2 

2 

3.08(-3) 

3.63C-2) 

°2 

32 

2.38(-2) 

1 . 76(-2) 

H20 

18 

2.14C-1) 

2.81(-1) 

N2 

28 

7.45(-l) 

6.28C-1) 

V.  Results  and  Discussion 


The  results  requested  in  the  B  problem  are  the  flame  velocity  and  the 
maximum  values  of  the  mass  fraction  profiles  of  H,  0,  OH,  and  HO2  and  the 
maximum  and  minimum  values  of  the  enthalpy.  The  mixture  enthalpy  profile  is 
computed  from  the  other  known  functions.  Specifically, 


mix 


N 

Z 

i=l 


Yihi 


(40) 


which  for  transport  method  5  with  constant  cp  reduces  to 
N 
=  I 


mix 


i=l 


Y.h.  +  c  (T-T  ). 
11  p  0 


(41) 


Since  the  functions  Y^(x)  and  T(x)  are  known  hm^x(x)  can  be  generated. 

The  flame  speeds  are  shown  in  Table  16.  The  first  column  contains  the 
values  predicted  by  the  code  using  the  nominal  values  of  i|>5,e  and  MINT. 

The  second  column  shows  the  values  predicted  when  the  burned  end  numerical 
boundary  iJir,  the  temporal  error  criterion,  e,  and  the  number  of  breakpoints 
NINT  are  simultaneously  changed.  The  third  column  shows  the  flame  speed  pre¬ 
dictions  from  the  same  code  and  input  parameters  except  that  the  kinetics  of 
Dixon-Lewis^  are  used. 

Since  our  grid^  adapts  to  the  temperature  profile,  the  flame  speed  from 
this  profile  is  the  most  accurate.  All  other  determinations  agree  within  5% 
except  the  HO2  profile.  The  reason  for  this  exception  is  that  the  interval  of 
integration  is  taken  from  the  cold  boundary  to  approximately  the  mid-point  of 
^i-space  (see  Eq.  (9)).  The  HO2  profile  (see  Fig.  1)  is  a  sharp  "spike"  that 
occurs  within  the  integration  interval  and  so  the  flame  speed  determination  is 
expected  to  be  grossly  in  error.  This  has  been  checked  by  raising  NCN  from  six 
to  twelve  and  simultaneously  lowering  ^5  to  2.0.  The  former  extends  the  region 
of  detail  about  the  inflection  of  the  temperature  profile  and  the  latter  de¬ 
creases  the  spatial  resolution.  In  this  case  the  flame  velocity  predicted  for 
the  HO2  profile  was  150  cm-s-1  showing  the  expected  trend.  In  addition  the 
flame  speed  predicted  for  the  OH  profile  was  197  cm-s”l  and  all  other  determi¬ 
nations  changed  by  3  units  or  less. 

All  these  results  are^in  good  agreement  with  the  computed  flame  velocity 
of  200  cm-s'l  (Dixon-Lewis’  his  figure  40).  The  corresponding  measurements’  show 
relatively  wide  scatter;  i.e.,  187-280  cm-s-1. 

Since  the  numerical  range  of  ip  space  is  finite,  we  need  to  check  the 
adequacy  of  our  selection  of  the  value  of  this  numerical  boundary.  If  this 
range  is  too  small,  mass  or  heat  or  both  will  "leak"  out. 

The  mass  and  heat  flux  fractions  are  defined 

Ek  =  Yk  +  pYk  Vk(pu)_1  (42) 
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TABLE  16.  COMPUTED  FLAME  SPEEDS 


Variable 

g 

Flame 

Velocity 

Flameb 

Velocity 

Flamec 

Velocit; 

cm-s"* 

cm-s~* 

-1 

cm-s 

Temp 

200 

200 

217 

H 

199 

200 

’6 

OH 

196 

190 

214 

0 

195 

197 

O 

H°2 

-1300 

-229 

121 

«2 

199 

200 

216 

C'J 

O 

199 

199 

216 

H20 

199 

199 

216 

N2 

197 

198 

216 

a*5  =  2.9,  e  =  3xl0"5, 

NINT  =  24, 

NCN 

=  6,  FC  =  6. 

b^5  =  3.3,  e  =  10"5,  NINT  =  30,  NCN  « 

6,  FC  =  6  (check) 

c  9 

Dixon-Lewis  Kinetics 

*5  =  0.8, 

e  = 

3xl0-4,  NINT  =  16, 

NCN  =  6,  FC  =  8, 

and 

e  *  Tu  ♦  A(T)x  cp_1  (pu)_1,  (43) 

respectively.  Since  Vr  is  always  related  to  the  derivative  (Yj()x  and  in  the 
Fick's  law  case,  considered  here,  (42)  can  be  written 

eK  *  YIC  -  »I’im<YK>x  •  <44> 


We  see  that  if  the  derivatives  (Yj^) x  *  (T)x  ®  0  at  the  cold  boundary  then  the 
flux  fraction  should  equal  the  originally  assigned  values  of  Yjg  and  Tu.  we 
find  the  derivatives  to  be  411  <_  10"?  and  that  the  values  of  the  flux  fractions 
are  identical  to  the  cold  boundary  values  to  at  least  five  significant  figures. 
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Table  17  shows  the  requested  peaks  for  the  Peters  B  problem.  These  are 
determined  from  the  output  routines  which  typically  contain  'v  200  data  points 
per  profile.  These  in  turn  have  been  generated  from  the  40  collocation  points 
and  the  known  coefficients  of  the  B-spline  functions. 


Mass  fraction  profiles  of  the  minor  (OH,  0,  H  and  HO2)  and  major  (H2,  O2, 
H2O  and  N2)  species  can  be  found  in  figures  1  and  2,  respectively.  Figure  3 
shows  the  profiles  of  temperature  and  enthalpy.  The  most  interesting  of  these 
profiles  are  those  of  HO2  and  OH.  The  fact  that  the  HO2  peak  lies  close  to  the 
smaller  OH  peak  suggests  that  they  may  be  related.  Let  us  first  examine  what 
is  happening  to  form  the  H02  peak. 

TABLE  17.  PEAK  VALUES  OF  RADICALS  AND  EXTREMA  OF  ENTHALPY 


i 

Y. 

1  max 

Location  (mm) 

Temp  (K) 

H 

2.674QE-3 

1.0935 

1353.8 

OH 

1.0012E-2 

2.4838 

1968.6 

0 

7 . 5460E-3 

1.1399 

1408.9 

ho2 

7.9263E-4 

0.78185 

369.0 

max  enthalpy 

8.7725E+1* 

0.90429 

916.9 

min  enthalpy 

-2.2690E+1* 

1.2143 

1479.5 

•units  of  cal/gm. 

The  only  source  of  H02  in  this  region  is  the  reaction 

H  ♦  02  ♦  M  ■*  H02  ♦  M  •  (R7) 

H  is  transported  toward  the  cold  boundary  by  diffusion  and,  before  the  flame 
front,  O2  is  in  abundant  supply.  So  H02  builds  up  to  the  HO2  peak  near  0.76  mm 
until  the  reactions 


H°2 

♦  H  -*•  2011 

(RS) 

ho2 

♦  H  -»  H2  ♦  02 

(R9) 

ho2 

♦  OH  H20  +  02 

(R10) 
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are  significant  and  together  dominate  the  HO2  loss.  (Small  but  non-negligible 
losses  to  the  HO2  concentration  also  occur  through  diffusion  and  conduction.) 
The  temperature  rises  sufficiently  till  the  rate  of  these  HO2  loss  mechanisms 
exceed  the  rate  of  HO2  production  and  Y^q  is  observed  to  fall. 

The  smaller  OH  peak  (near  0.8  mm)  is  caused  by  the  following  mechanism. 

The  major  source  of  OH  is  reaction  8  (which  also  depletes  H02);  the  major 
OH  sink  is 


H2  +  OH  H20  +  H, 


(Rl) 


but  the  OH  profile  tends  to  track  the  HO.,  profile  except  that  near  0.9  mm  the 
chain  branching  reactions 


02  +  H  +  0  +  OH 
H2  +  0  ■>  H  +  OH 


(R2) 

(R3) 


become  important.  With  increasing  temperature  the  OH  builds  up  again.  The 
gap  in  the  Yqh  profile  is  the  location  where  OH  is  tracking  H02  and  reactions 
(R2)  and  (R3)  do  not  have  significant  influence. 

In  Figure  4  we  see  a  plot  of  Xn^  vs.  distance.  The  rise  needs  to  be 
understood  because  N2  is  a  non-reacting  diluent  and  the  corresponding  Yjj  plot 
is  flat  (see  Fig.  2).  As  we  shall  now  show,  this  behavior  can  be  accounted 
for  in  the  Y  to  X  transformation. 


Given  that  Yi(A)  =  Yj (B)  =  ...,  where  A,B, ...  are  spatial  locations,  we 
hat  is  the  relationship  between  say,  Xi(A)  and  X^B).  In  general 

X^A)  =  (Y.(A)/M.)/(E  Y.(A)/M.).  (45) 


Then 


Xi (B)  =  (Y^BJ/M.J/CE  Yj  (B)/Mj ) 


(46) 


=  (Yi(A)/M.)/(Z  Y. (B)/M.)  . 

We  require  an  expression  for  the  denominator  of  (46)  in  terms  of  location  A. 
Since 

£  Y.(B)/M  =  E  Y.(B)/M  +  Y  (B)/M 

J  J  j^i  1  )  1  1 

and 


E  Y. (A)/M.  =  E  Y.(A)/M.  +  Y. (A)/M.  . 
J  1  1  J  11 
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HOLE  FRAC 


Fig.  4  N„  mole  fraction  profile. 


We  can  equate  Y^(B)  and  Y^(A)  to  find  that 


I  Y j  (B) /M^  =  I  (Yj  (B)  -Yj  (Ajj/Mj  +  Z  Y ^  (A) /M^  (47) 

Substituting  (47)  into  (44)  and  using  (45)  we  find  that 


X.(A)/X.(B)  -  1  +  a 


(48) 


where  a  =  Z  ( (Y .(B) -Y  (A) )/M  )/ (Z  Y  (A)/M  ). 

j?fi  3  J  J  ■>  > 

Note  that  if  a  <<  1  then  X. (A)  =  X.(B);  if  a  >  0  then  X. (A)  >  X.(B).  These  last 
inequalities  allow  the  possibility  that  a  flat  Y.  profile  couli  be  transformed 
into  an  X^  profile  that  is  quite  different  in  shape. 

As  a  check,  we  take  the  Yj 's  from  the  code  output  and  compute  e  (see  Eq. 
44),  for  a  selected  number  of  collocation  points  16  through  21  corresponding 
to  0.7  mm  to  0.8  mm.  We  use  a  point  in  the  flat  portion  of  the  curve,  location 
10  (s*  0.52  mm)  as  the  fiducial  or  reference  point,  A.  From  the  code  output  we 
can  also  compute  X.(A)/X(B)  -  1.  If  nothing  but  the  transformation  were  in¬ 
volved  these  values  woulcl  be  identical.  Table  18  shows  these  values.  As  can 
be  seen  the  values  are  quite  close;  the  differences  are  attributed  to  transport 
processes  and  numerical  "slop".  ^ 

(Recall  that  Y^  is  computed  from  1-E  Y^.) 


Figure  5  shows  Y0, 


as  a  function  of  distance  and  the  expected  monatomic 


decrease.  However  Figure  6  shows  Xq2  and  a  sensible  "bump"  appears  right  before 
the  steep  slope.  Reference  to  Table  18  shows  that  at  locations  16  5  17  the 
trend  is  similar  to  the  Yj^.  However  for  the  higher  locations  the  actual  curve, 
G,  departs  from  the  one  that  only  takes  into  account  the  Y-*X  transformation,  F. 
In  other  words  some  other  process  is  occurring  besides  the  simple  transformation 
Examination  of  the  output  shows  that  O2  diffusion  is  the  major  influence  on 
the  profile.  That  reactions  play  a  negligible  part  in  explaining  the  Yj^  "bump" 
is  consistent  with  the  low  temperatures  in  this  region  (<500K) . 


VI.  SUMMARY 


Besides  the  requested  information  we  have  found  several  interesting  profiles 
whose  bumps  and  wiggles  could  have  been  attributed  to  numerical  artifacts,  but 
upon  closer  examination  proved  to  be  real  for  this  flame  and  its  input  parameters 

In  addition,  we  provide  a  set  of  nonkinetic  input  parameters  for  this 
stoichiometeric  hydrogen-air  flame  and  compare  these  values  against  known 
experimental  determinations.  Finally,  we  show  in  the  Appendix  that  the  multi¬ 
species  transport  algorithm  is  not  important  unless  one  is  considering  radical 
profile  accuracies  of  s/20%  or  less.  We  also  show  that  the  major  differences 
between  the  two  algorithms  employed  can  be  attributed  to  the  effects  of  thermal 
diffusion. 
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TABLE  j 

L8.  VALUES  OF 

Y-*X  TRANSFORMATION  AT 

SPECIFIC  LOCATIONS 

Location 

(-B> 

F(N2,B,10) 

G(N2>B,10) 

F(02,B,10) 

G(02,B,10) 

16 

-1.80(-3) 

-1.90(-3) 

-1.66(-3) 

-1.62 (-3) 

17 

-4.10(-3) 

-4. 30 (-3) 

-3. 77 (-3) 

-3.30(-3) 

18 

-6.99(-3) 

-7.79(-3) 

-6.36(-3) 

-4. 78 (-3) 

19 

-1.28(-2) 

-1.31 (-2) 

-1.12(-2) 

-3.91 (-3) 

20 

-1.83(-2) 

-1.85(-2) 

-1.14(-2) 

♦3.73(-3) 

21 

F(i,B,A)  =  S 

G(i,B,A)  * 

-2.68(-2)  -2.69C-2) 

((Y. (B)  -  Yj(A)/Mj)/(E  Y^(A)/M^) 

CX.(A)/Xi(B))  -  1 

-1.84(-2) 

+3.05C-3) 
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APPENDIX  A 


Comparison  of  Results  Using  Transport  Method  6. 

In  the  main  body  of  the  text  we  have  employed  a  mixing  algorithm  which  is 
based  upon  Equations  (29),  (30)  and  (39).  This  method,  our  so-called  method 
5, 20  consists  of  a  simple  definition  for  Xmix,  Eq.  (29),  and  a  generalization 
of  Pick's  Law  plus  the  selection  of  constants  for  p2  Dim,  p\  and  c„.  In  order 
to  check  the  effects  of  this  rather  simple  mixing  algorithm  we  employed  our 
so-called  method  6.  This  method  also  used  Eq.  (29)  as  the  definition  of  Xmix 
but  accounts  for  the  diffusion  velocity  by  solving  a  modified  form  of  the 
Stefan-Maxwell  equations.  In  this  method  the  transport  coefficients  are 
temperature  dependent.  This  method  also  crudely  accounts  for  thermal  diffusion. 
We  have  assumed  that  only  H  and  H2  have  non-zero  thermal  diffusion  coefficients. 
Futher  details  are  given  in  Ref.  20.  All  other  parameters  remain  the  same. 

Table  A-l  shows  the  computed  flame  speed  for  both  the  "Peters  B"  kinetic 
listed  in  Table  1  and  the  Dixon-Lewis  kinetics  listed  in  Table  2.  These  values 
are  about  4%  greater  than  those  in  Table  16.  The  exceptions  are  for  HO2  whose 
difficulties  we  have  discussed  above  and  for  N2.  The  Yn2  are  practically 
straight  lines  of  constant  value  and  so  the  use  of  Eq.  (9)  to  obtain  flame 
speeds  will  not  in  general  be  reliable. 

Figures  A-l  and  A-2  show  the  minor  (H,0H,0  and  HO2)  and  the  major  (H2,02 
H2O  and  N2)  species  profiles  respectively.  Figure  A-3  shows  the  temperature 
and  enthalpy  profiles  for  transport  method  6.  Qualitatively  there  is  little 
difference  between  these  profiles  and  the  corresponding  ones  of  method  5, 

Figures  1,2,  and  3  of  the  text.  Figure  A-4  shows  three  plots  of  the  enthalpy, 
the  function  with  the  greatest  observed  difference.  The  solid  line  is  the 
enthalpy  computed  using  transport  method  5,  the  dotted  lines  is  method  6  and 
the  dash-dot  line  is  method  6  with  the  thermal  diffusion  set  equal  to  zero. 

As  can  be  seen  almost  all  the  difference  between  method  5  and  method  6  can  be 
attributed  to  the  effects  of  thermal  diffusion. 

T 

We  add  in  passing  that  for  transport  method  6  with  D  =  0  the  temperature 
and  major  species  O2,  H2O  and  N2  practically  overlaid  the  corresponding  plots 

of  both  method  5  and  method  6  (D‘  f  0) .  For  H-  method  6  showed  YH2  to  decline 

at  smaller  values  of  x  while  methods  5  and  6  (uT  =  0)  nearly  overlapped. 

The  minor  species  profiles  differences  are  characterized  by  the  Yi  peak 
values  as  shown  in  Table  A-2.  At  the  peak  and  in  the  post -flame  region  the 
values  for  the  minor  species  computed  using  method  6(D^=0)  were  larger  than  the 
values  computed  using  method  6(DT  /  0)  and  these  in  turn  were  larger  than  the 
values  computed  using  method  5. 

T 

The  flame  speed  computed  from  the  temperature  profile  for  method  6(D  =  0) 

was  230  cms  .  This  is  about  12%  greater  than  the  results  of  method  6(DT  j<  0). 

Wamatz's22  computations  are  for  a  45%  H2  in  air  premixed  flame.  In  this  case  he 
finds  that  setting  =  0  raises  the  flame  speed  by  'v  5%. 
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Wamata,  J.  Ber.  Bun8engee.  Phys.  Chem.  82,  643-649, 


(1978). 


i 


TABLE  A-l .  COMPUTED  FLAME  SPEEDS  METHOD  VI 


Variable 

Flamea 

Velocity 

Flame*5 

Velocity 

Temp. 

208  cm/s 

215  cm/s 

H 

207 

221 

OH 

193 

222 

0 

203 

221 

ho2 

-150 

-3638 

H2 

207 

215 

°2 

207 

215 

h2o 

207 

215 

N2 

205 

154 

a)  <l>  =  2.9,  e  =  3xl0-5,  NINT  =  24,  NCN  =  6,  FC  =  6. 

b)  Dixon-Lewis  Kinetics9:  <J»g  =  2.9,  e  *  3x10  5,  NINT  «  24,  NCN  =  6,  FC 
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TABLE  A-2 .  COMPARISON  OF  RADICAL  PEAK  HEIGHTS  AND  ENTHALPY 

TWO  TRANSPORT  ALGORITHMS. 


Species 

Function 

Method  5 

Method  6 

H 

Yi(max) 
Temp.(K) 
location (mm) 

2.6740E-3 

1353.8 

1.0935 

2.6852E-3 

1369.3 

1.1051 

OH 

Y 

i (max) 
Temp(K) 
location (mm) 

1.0012E-2 

1968.6 

2.4838 

1.0356E-2 

1973.0 

2.4878 

0 

Yi(max) 

Temp. (K) 
location (mm) 

7.5460E-3 

1408.9 

1.1399 

7.8953E-3 

1427.7 

1.1567 

H02 

Yi (max) 

Temp. (K) 
location (mm) 

7.9263E-4 

369.0 

0.78185 

6.6590E-4 

382.2 

0.76351 

enthalpy* 

maximum 

Temp. (K) 
location (mm) 

8.7725E+1 

916.9 

0.90429 

6.6526E+1 

901.8 

0.89672 

enthalpy* 

minimum 

Temp. (K) 
location (mm) 

-2.2690E+1 

1479.5 

1.2143 

-3.2441E+1 

1469.3 

1.2005 

*  units  of  cal/gm 


EXTREMA  FOR 


Method  6(DT=0) 


2.8494E-3 

1364.6 
1.1049 

1.0409E-2 

1965.8 

2.5436 

7.9758E-3 

1499.1 

1.2452 

7.0769E-4 

377.2 

0.75791 

8. 1928E+1 
928.5 
0.90314 

-2.3675E+1 

1484.7 
1.2269 
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Fig.  A-2. 


0*0091  0*0001  0*000  0*0  0*000-  0*0001 

(H0/TbO)H  HO  (X)dW31 
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fig.  A-4.  Enthalpy  profits;  tlethod  5  M_ 


GLOSSARY 


=  heat  capacity  at  constant  pressure,  cal-mole~*-K”* 

c  =  c  .Y.,  the  specific  heat  of  the  mixture,  cal-em’1-K"^ 
p  pi  1  r  ° 

Cp^  =  specific  heat  of  ith  species,  cal-gm"*-K~* 

2  -1 

=  binary  diffusion  coefficients  for  species  k  and  j ,  cm  -s 
lip0  =  enthalpy,  cal-mole*^ 
hp  =  specific  enthalpy,  cal-gm”^ 
k  =  Boltzmann's  constant,  =  1.3306  erg  -  K~* 


k.  =  rate  coefficient  for  the  ith  reaction,  in  centimeter-mole-seconds 

l  * 

units 


N  =  number  of  species 

Nh  =  molecular  weight  of  ith  species,  gm-mole~* 

n  -2  -1 

mQ  =  mass  flux,  gm-cm  -s 
P  =  total  pressure,  atmosphere 

R  =  Gas  constant,  =  1. 9372  cal-mole  ^K-1,  =  82.05  cm3-atmos-mole~1'-K'1 

_3 

Rp  =  Rate  of  production  of  ith  species  by  chemical  reactions,  mole-cm  -s 
=  burning  velocity,  cm-s 
T  =  temperature, K 


Tg  =  temperature  of  burned  mixture,  K. 

T*  =  TCc./K)'1 

Tu  =  temperature  of  unburned  mixture,  K 
t  =  temporal  coordinate,  s 

u  =  fluid  velocity,  cm-s-1 

Vp  =  diffusion  velocity  of  ith  species,  cm-s 


Y. 

l 


=  mole  fraction  of  ith  species 
=  mass  fraction  of  ith  species 
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GLOSSARY  'cont.) 


Y^B  =  mass  fraction  of  ith  species  in  the  burned  mixture. 

Y^u  =  mass  fraction  of  ith  species  in  the  unbumed  mixture. 

e^/k=  L-J  or  Stockmayer  parameter,  K 

=  viscosity  of  ith  ;pecies,  gm-s”*-cm* 

X  =  mixture  conductivity,  cal-cm’^-s'^-K-* 

X^  =  conductivity  of  ith  species,  cal-cm~*-s”*-K-* 

-8  ° 

=  L-J  or  Stockmayer  parameter,  10”  cm  (=  A) 

_3 

p  =  fluid  density,  gm-cm 

_2 

ij;  =  transformation  coordinate  (see  equation  5) ,  gm-cm 
12  2  2 

n  ’  ,0  *  =  collision  integrals  (tabulated  in  reference  1) 
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Monterey,  CA  93940 
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2  Commander 

Naval  Ordnance  Station 
ATTN:  S.  Mitchell 
Tech  Lib 

Indian  Head,  MD  20640 

2  AFOSR 

ATTN:  L.  Caveny 

B.T.  Wolfson 
Bolling  AFB,  DC  20332 

2  AFRPL  (DYSC) 

ATTN:  D.  George 

J.N.  Levine 

Edwards  AFB,  CA  93523 

2  National  Bureau  of  Standards 
ATTN:  J.  Hastie 

T.  Kashiwagi 
Washington,  DC  20234 

1  Lockheed  Missile  6  Space  Co. 
ATTN:  Tech  Info  Ctr 
3521  Hanover  Street 
Palo  Alto,  CA  94304 

1  Aerojet  Solid  Propulsion  Co. 
ATTN:  P.  Micheli 
Sacramento,  CA  95813 

1  ARO  Incorporated 
ATTN:  D.  Dougherty 
Arnold  AFB,  TN  37389 

1  Atlantic  Research  Corporation 
ATTN:  M.K.  King 
5390  Cherokee  Avenue 
Alexandria,  VA  22314 


2  Cal span  Corporation 
ATTN:  E.B.  Fisher 
A.P.  Trippe 
P.0.  Box  400 
Buffalo,  NY  14225 

1  01 in  Corporation 

Badger  Army  Ammunition  Plant 
ATTN :  J .  Ramnarace 
Baraboo,  WI  53913 

1  Foster  Miller  Associates 
ATTN:  A.J.  Erickson 
135  Second  Avenue 
Waltham,  MA  02154 

1  General  Electric  Company 
Armament  Department 
ATTN:  M.J.  Bulman 
Lakeside  Avenue 
Burlington,  VT  05402 

1  General  Electric  Company 
Flight  Propulsion  Division 
ATTN:  Tech  Lib 
Cincinnati,  OH  45215 

1  General  Motors  Research 

Laboratories 
Physics  Department 
ATTN :  C .  Dasch 
Warren ,  MI  48090 

2  Hercules  Powder  Co. 

Allegheny  Ballistics  Lab. 
ATTN:  R.  Miller 

Tech  Lib 

Cumberland,  MD  21501 


1  AVCO  Everett  Research  Lab 
ATTN:  D.  Stickler 
2385  Revere  Beach  Parkway 
Everett,  MA  02149 


1  Hercules  Incorporated 
Bacchus  Works 
ATTN:  B.  Isom 
Magna,  UT  84044 
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I  Paul  Gough  Associates,  Inc. 
ATTN:  P.  S.  Gough 
P.  0.  Box  1614 
Portsmouth,  NH  03801 

1  Physics  International  Company 
2700  Merced  Street 
Leandro,  CA  94577 


1  Pulsepower  Systems,  Inc. 
ATTN:  L.  C.  Elmore 
815  American  Street 
San  Carlos,  CA  94070 

3  Rockwell  International  Corp. 
Rocketdyne  Division 
ATTN :  C .  Obert 

J.  E.  Flanagan 
A.  Axeworthy 
6633  Canoga  Avenue 
Canoga  Park,  CA  91304 


1  Hercules,  Inc. 

Industrial  Systems  Dept 
P.  0.  Box  548 
ATTN:  W.  Haymes 
Tech  Lib 

McGregor,  TX  76657 
1  SAI 

ATTN:  R.  Farmer 
21133  Victory  Boulevard 
Suite  216 

Canoga  Park,  CA  91303 

1  Science  Applications,  Inc. 
ATTN:  R.  B.  Edelman 
Combustion  Dynamics  5 
Propulsion  Division 
23146  Cumorah  Crest 
Woodland  Hills,  CA  91364 


No.  of 

Copies  Organization 

1  Shock  Hydrodynamics,  Inc. 
ATTN:  W.  H.  Anderson 
4710-16  Vineland  Avenue 
North  Hollywood,  CA  91602 

1  Thiokol  Corporation 
Elkton  Division 
ATTN:  E.  Sutton 
Elkton,  MD  21921 

3  Thiokol  Corporation 
Huntsville  Division 
ATTN:  D.  Flanigan 
R.  Glick 
Tech  Lib 

Huntsville,  AL  35807 

2  Thiokol  Corporation 
Wasatch  Division 
ATTN:  J.  Peterson 

Tech  Lib 
P.  0.  Box  524 
Brigham  City,  UT  84302 


1  TRW  Systems  Group 
ATTN:  H.  Korman 
One  Space  Park 
Redondo  Beach,  CA  90278 

2  United  Technologies 
Chemical  Systems  Division 
ATTN:  R.  Brown 

Tech  Lit 

Sunnyvale,  CA  94086 

1  Universal  Propulsion  Co. 

ATTN:  H.  J.  McSpadden 
Black  Canyon  Stage  1,  Box  1140 
Phoenix,  AZ  85029 
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1  IITRI 

ATTN:  M.  J.  Klein 
10  West  35th  Street 
Chicago,  IL  60616 

1  Battelle  Memorial  Institute 
ATTN :  Tech 
505  King  Avenue 
Columbus,  OH  43201 


1  Brigham  Young  University 

Dept  of  Chemical  Engineering 
ATTN:  M.  W.  Beckstead 
Provo,  UT  84601 


1  California  Institute  of  Tech 
204  Kamar  Lab 
Mail  Stop  301-46 
ATTN:  F.  E.  C.  Culick 
1201  E.  California  Street 
Pasadena,  CA  91125 


1  Case  Western  Reserve  Univ- 

Division  of  Aerospace  Sciences 
ATTN:  J.  Tien 
Cleveland,  OH  44135 

3  Georgia  Institute  of  Technology 
School  of  Aerospace  Engineering 
ATTN:  B.  T.  Zinn 
E.  Price 
W.  C.  Strahle 
Atlanta,  GA  30332 

1  Institute  of  Gas  Technology 
ATTN:  D.  Gidaspow 
3424  S.  State  Street 
Chicago,  IL  60616 


No.  of 

Copies  Organization 

1  Johns  Hopkins  University/APL 

Chemical  Propulsion  Information 
Agency 

ATTN:  T.  Christian 
Johns  Hopkins  Road 
Laurel,  MD  20707 

1  Massachusetts  Institute  of 

Technology 

Dept  of  Mechanical  Engineering 
ATTN:  T.  Toong 
Cambridge,  MA  02139 

1  Princeton  Combustion  Rsch  Center 

ATTN:  M.  Summerfield 
1041  US  Highway  One  North 
Princeton,  NJ  08540 

1  Pennsylvania  State  University 

Applied  Research  Lab 
ATTN:  G.  M.  Faeth 
P.  0.  Box  30 

State  College,  PA  16801 


1  Pennsylvania  State  University 

Dept  of  Mechanical  Engineering 
ATTN:  K.  Kuo 

University  Park,  PA  16802 

1  Pennsylvania  State  University 

Dept  of  Material  Sciences 
ATTN:  H.  Palmer 
University  Park,  PA  16802 

1  Princeton  University 

Forrestal  Campus 
ATTN:  I.  Glassman 
Tech  Lib 
P.  0.  Box  710 
Princeton,  NJ  08540 
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2  Purdue  University 

School  of  Mechanical  Engineering 
ATTN :  J .  Osborn 

S.  N.  B.  Murphy 
TSPC  Chaffee  Hall 
West  Lafayette,  IN  47906 

1  Rutgers  State  University 
Dept  of  Mechanical  and 
Aerospace  Engineering 
ATTN:  S.  Temkin 
University  Heights  Campus 
New  Brunswick,  NJ  08903 

1  Southwest  Research  Institute 
Fire  Research  Section 
ATTN:  W.  H.  McLain 
8500  Culebra  Road 
San  Antonio,  TX  78228 

4  SRI  International 
ATTN:  Tech  Lib 
J.  Barker 
D.  Corsley 
D.  Golden 

333  Ravenswood  Avenue 
Menlo  Park,  CA  94025 


1  Stevens  Institute  of  Technology 
Davidson  Laboratory 
ATTN:  R.  McAlevy,  III 
Hoboken,  NJ  07030 


1  University  of  California, 
San  Diego 
Ames  Department 
ATTN:  F.  Williams 
P.  0.  Box  109 
La  Jolla,  CA  92037 


1  University  of  Illinois 

Dept,  of  Mechanical  Engineering 
144  MEB,  1206  W.  Green  St. 

Urbana,  IL  61801 

1  University  of  Minnesota 

Dept  of  Mechanical  Engineering 
ATTN:  E.  Fletcher 
Minneapolis,  MN  55455 

1  University  of  Southern  California 

Dept  of  Chemistry 
ATTN:  S.  Benson 
Los  Angeles,  CA  90007 


2  University  of  Texas 

Dept  of  Chemistry 
ATTN:  W.  Gardiner 
H.  Schaefer 
Austin,  TX  78712 


1  Rensselaer  Polytechnic  Institute 

ATTN:  Prof.  A.  Fontijn 
Troy,  N.  Y.  12181 

3  Center  for  Non  Linear  Studies 

ATTN:  Dr.  Basil  Nichols 
Mail  Stop  457 

Los  Alamos  National  Laboratory 
Los  Alamos,  New  Mexico  87545 

1  Princeton  University 

Department  of  Mechanical  and 
Aerospace  Engineering 
ATTN:  Prof.  Rolf  D.  Reitz 
Princeton,  NJ  08540 


2  University  of  Utah 

Dept  of  Chemical  Engineering 
ATTN:  A.  Baer 

G.  Flandro 

Salt  Lake  City,  UT  84112 
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Aberdeen  Proving  Ground 

Dir,  USAMSAA 
ATTN:  DRXSY-D 

DRXSY-MP,  H.  Cohen 

Cdr,  USATECOM 

ATTN:  DRSTE-TO-F 

Dir,  USCSL,  Bldg.  E3516,  EA 
ATTN:  DRDAR-CLB-PA 
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